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Abstract 


This  paper  describes  the  software  for  a  new  modified  Cholesky  factorization  recently  proposed  by  the 
authors.  Given  a  symmetric  but  not  necessarily  positive  definite  matrix  A,  the  modified  Cholesky  factor¬ 
ization  computes  a  Cholesky  factorization  of  A  4-  E ,  where  E  —  0  if  A  is  safely  positive  definite,  and  E  is 
a  diagonal  matrix  chosen  to  make  A  -r  E  positive  definite  otherwise.  The  modified  Cholesky  factorization 
was  introduced  by  Gill  and  Murray  and  refined  by  Gill,  Murray  and  Wright,  and  is  commonly  used  in 
optimization  algorithms.  Our  version,  which  is  based  upon  new  techniques,  has  a  considerably  smaller 
a  priori  upper  bound  on  the  size  of  E  than  the  Gill.  Murray  and  Wright  factorization,  and  appears  to 
generally  produce  a  smaller  £,  and  a  well-conditioned  A  +  E,  in  practice.  Its  cost,  like  the  Gill,  Murray 
and  Wright  version,  is  only  a  small  multiple  of  n2  operations  greater  than  the  standard  Cholesky  factor¬ 
ization.  Thus  it  may  be  useful  in  optimization  algorithms.  We  summarize  o"r  algorithm  and  describe 
the  code  and  its  use. 


1  Introduction 


This  paper  describes  the  software  for  a  new  modified  Cholesky  factorization  algorithm  presented  in 
Schnabel  and  Eskow  [19SS] -  The  modified  Cholesky  factorization  is  intended  for  matrices  that  are  sym¬ 
metric  but  not  necessarily  positive  definite.  Given  a  matrix  A  6  Rnxn.  the  modified  Cholesky  factoriza¬ 
tion  computes 

Pt(A  +  E)P  =  LLT(orLDLT). 

where  P  is  a  permutation  matrix,  and  E  €  Rnxn  is  0  if  .4  is  safely  positive  definite,  otherwise  £  is  a 
nonnegative  diagonal  matrix  chosen  such  that  AA-  E  is  safely  positive  definite.  This  type  of  factorization 
was  introduced  by  Gill  and  Murray  [1974]  and  subsequently  refined  by  Gill,  Murray,  and  Wright  [19!?  1]. 
It  has  become  important  in  solving  problems  in  optimization,  and  is  used  in  many  line  search  methods 
for  unconstrained  and  constrained  optimization  problems  (Gill.  Murray,  and  Wright  [1981])  as  well  as  in 
some  trust  region  methods  (Dennis  and  Schnabel  [1983]). 

The  modified  Cholesky  factorization  of  Schnabel  and  Eskow  has  superior  theoretical  properties 
to  the  method  of  Gill.  Murray,  and  Wright,  and  appears  to  have  computational  advantages  as  well.  The 
upper  bound  on  j|£||M  for  Schnabel  and  Eskow's  method  is  at  worst  about  2na,  where  a  =  max  S.4,y  |.  as 
opposed  to  roughly  n- a  for  the  method  of  Gill,  Murray,  and  Wright.  In  addit  ion,  extensive  computational 
testing  of  the  software  described  in  this  paper  was  performed  by  Schnabel  and  Eskow  [1988],  on  a  set 
of  randomly  generated  indefinite  matrices  with  n  =  25.  50  and  75.  The  norm  of  the  matrix  £  and 
condition  number  of  the  (.4  -t-  £)  computed  by  this  algorithm  were  compared  to  those  produced  by 
the  Gill.  Murray,  and  Wright  [1981]  factorization.  In  almost  all  cases,  ||£||^  for  the  new  algorithm 
was  significantly  smaller,  while  both  methods  consistently  produced  acceptably  conditioned  matrices. In 
addition,  the  ||£||(x  produced  by  the  Schnabel  and  Eskow  algorithm  almost  always  was  within  a  factor  of 
2  of  the  magnitude  of  the  most  negative  eigenvalue  of  A,  and  often  much  closer.  Another  nice  property 
of  the  new  factorization  (Gould  [1989])  is  that  in  applications  like  multifrontal  methods,  where  .4  is  large 
and  sparse  and  no  pivoting  is  performed,  the  matrix  A  may  be  generated  and  processed  incrementally, 
as  the  factorization  proceeds,  whereas  the  Gill.  Murray,  and  Wright  method  requires  the  entire  matrix 
to  be  known,  and  processed,  in  its  initialization  stage.  Thus,  the  new  factorization  may  be  useful  m 
optimization  and  other  contexts. 

In  Section  2,  we  give  a  brief  description  of  the  algorithm.  Section  3  demonstrates  the  use  of  the 
algorithm  on  a  small  (4  x  4)  example  matrix  and  section  4  explains  the  parameters  and  organization  of 
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the  code  Appendices  A  and  B  provide  a  sample  driver  and  its  output,  respectively,  and  Appendix  t' 
contains  the  .ode  for  the  modified  Cholesky  factorisation. 

2  The  Modified  Cholesky  Factorization  Algorithm 

This  section  briefly  describes  the  modified  Cholesky  factorization  algorithm  presented  in  Schnabel  and 
Eskow[198$E  Some  detail  concerning  the  techniques  used  to  prevent  an  ill-conditioned  result  is  included 
because  the  user  has  the  option  of  adjusting  two  tolerances  related  to  this.  The  cost  of  the  factorization 
is  also  described.  For  a  more  detailed  explanation  of  the  entire  algorithm,  see  the  original  paper. 

The  new  modified  Cholesky  factorization  uses  a  two  phase  approach  to  compute  Pr(A  —  E)P  — 
LLJ .  where  .4  €  Rnxn  is  symmetric.  E  is  a  nonnegative  diagonal  ma'rix  or  0.  L  is  a  lower  triangular 
matrix  and  P  is  a  permutation  matrix.  Phase  1  computes  the  standard  Cholesky  factorization,  and 
ends  when  the  next  iteration  of  the  standard  factorization  would  cause  some  diagonal  element  in  the 
remaining  submatrix  to  become  non-positive.  It  performs  diagonal  pivoting  based  on  the  maximum 
diagonal  element  of  the  submatrix  remaining  to  be  factored.  Schnabel  and  Eskow  show  that  at  the  end 
of  this  phase,  no  element  in  the  submatrix  that  remains  to  be  factored  is  larger,  in  magnitude,  than  the 
sum  of  the  magnitudes  of  the  largest  diagonal  and  off-diagonal  elements  in  that  submatrix  originally. 

Phase  2  does  a  modified  Cholesky  factorization,  meaning  that  at  each  iteration,  an  amount 
Ejj  >  0  is  added  to  the  pivot  A,}  before  the  elimination  step.  A  diagonal  pivoting  strategy  is  also  used, 
but  in  this  phase  the  pivot  row  is  chosen  to  be  the  one  with  the  maximum  lower  Gerschgorin  bound 
estimate  (see  below).  The  amount  E;;  added  to  the  pivot  element  at  each  iteration  is  determined  from 
the  actual  lower  Gerschgorin  bound  of  the  pivot  row.  Schnabel  and  Eskow  show  that  due  to  the  choice 
of  E}j.  the  Gerschgorin  bounds  of  the  next  remaining  submatrix  (after  elimination)  do  not  grow.  In 
turn,  this  implies  that  || £!!<*,  is  bounded  above  by  the  magnitude  of  the  most  negative  lower  Gerschgorin 
bound  of  the  submatrix  that  remained  to  be  factored  at  the  start  of  phase  2.  This,  combined  with  the 
growth  bound  for  phase  1  leads  to  the  theoretical  result  mentioned  in  Section  1. 

The  entire  algorithm  is  outlined  in  Algorithm  2.1  below. 
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Algorithm  2.1  -  Modified  Cholesky  Decomposition 


Given  A  £  Jlnxn  symmetric  and  rl  and  r2  (e.g.  =  r2  =  macheps  j). 

find  factorization  LLT  of  .4  +  E,  E  >  0 


:=  max  ld.il 

1  <*<n 


J  :=  1 


(*  Phase  One,  A  potentially  positive  definite  *) 

While  j  <  n  do 

Pivot  on  maximum  diagonal  of  remaining  submatrix 


If 


min  {A„  -  -ii}  <  rp 

j  +  1  <i<n 


then  go  to  Phase  Two 

else  perform  j<|l>  iteration  of  standard  Cholesky  factorization 
and  increment  j 

(*  Phase  Two.  A  not  positive  definite  *) 

k  :=  j  —  1  (*  k  =  number  of  iterations  performed  in  Phase  One  *) 

Calculate  lower  Gerschgorin  bounds  of  d*  +  i 
For  j  :=  k  +  1  to  n  —  2  do 

Pivot  on  maximum  lower  Gerschgorin  bound  estimate 
Calculate  and  add  Ej3  to  dj; 


(*  Ejj  =  max{0.-d;;  +  max{  I'M.  ’V)'}.  £>-i.;-i}  *) 

•  =j  +  i 

update  Gerschgorin  bound  estimates 
perform  j‘h  iteration  of  Cholesky  factorization 
complete  factorization  of  final  2x2  submatrix  using  its  eigenvalues 


In  order  to  prevent  A  -h  E  from  being  singular  or  very  ill-conditioned,  the  algorithm  includes  the 
following  details.  Let  y  be  the  maximum  diagonal  element  of  A,  and  Tj  and  r2  be  some  small  constants 
(our  default  choice  is  r.  =  r2  =  macheps j).  The  switch  to  phase  2  is  made  when  some  diagonal  element 
in  the  remaining  submatrix  would  become  less  than  ny,  implying  that  only  an  indefinite  matrix  or  a 
positive  definite  matrix  with  condition  number  greater  than  T.  may  be  perturbed.  During  phase  2,  the 
the  amount  E}j  to  add  to  Ajj  is 

n 

=  max{0, -Ajj  +  max{  ^  \Atj |,  r2y}, £y_i,>_i} 

*=;  +  i 

where  £y_iy_i  is  the  amount  added  to  Aj-ij-i  in  the  previous  iteration.  The  r2y  term  in  the  above 
computation  allows  the  condition  number  of  A  +  E  to  be  bounded  above.  The  final  place  in  the  algorithm 
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where  the  conditioning  of  the  resultant  matrix  is  addressed  is  in  the  final  ((n  —  1)J|)  iteration  of  the 
factorization.  Eigenvalues  X!u  and  of  the  remaining  2x2  submatrix  are  computed,  which  are  th<>n 
used  to  calculate 

=  En  n  =  max{0,  En_;i0_2.  -A;o  +  r2  *  max{  - -j- -(Afc,  -  ), ->-}}. 

This  causes  the  E  norm  of  the  resultant  final  2x2  submatrix  to  be  no  greater  than  and  in  practice 
usually  results  in  En_i  n_i  ha"ing  a  smaller  value  than  would  otherwise  be  obtained  using  Gerschgorin 
bounds.  The  analysis  in  Schnabel  and  Eskow  [1988]  includes  these  details.  In  practice  the  condition 
number  of  A  4-  E  is  usually  no  greater  than  10/  minlrt.ro},  although  this  bound  does  not  hold  in  theory 
Phase  2  of  the  algorithm  pivots  on  an  estimate  of  the  lower  Gerschgorin  bounds  of  the  remaining 
submatrix.  Let  G,.{j  <  i  <  n)  denote  the  lower  Gerschgorin  bound  estimates  used  during  the  y 
iteration.  The  actual  lower  Gerschgorin  Circle  Theorem  bounds  are  computed  once  at  t he  start  of  phase 
2.  giving 

i  - 1  n 

G,  =  -d„  -  y:  |-4, •  k  |  -  y  M*i|,  i  =  j.  -  ■  ■ .  n. 

k=j  t=>  +  l 

where  j  o  the  iteration  in  which  the  algorithm  switches  to  phase  2.  Thereafter,  at  eacli  iteration  j  the 
bounds  are  estimated  by 

M,;|  Mo- I 

c,  =0.4- Mol - ^ - ,  i  =  n. . 

AJJ 

n 

Since  in  the  calculation  of  the  estimates  of  the  bounds,  the  sum  ^  Mol  needs  only  to  be  computed 

i=;  +  l 

once  at  each  iteration,  the  cost  of  computing  the  bound  estimates  is  at  most  n- j 2  each  additions  and 
multiplications  over  the  entire  algorithm,  whereas  it  would  be  0(n3)  if  the  actual  Gerschgorin  bounds 
were  used.  The  cost  of  the  entire  modified  Cholesky  factorization  is  at  most  2n2  additions  and  rr / 2 
multiplications  greater  than  the  n3/6-t-0(n)  each  multiplications  and  additions  for  the  standard  Cholesky 
factorization  of  positive  definite  matrices.  If  A  is  safely  positive  definite,  there  is  no  extra  cost. 
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3  Example  using  the  Modified  Cholesky  Factorization 


The  following  discussion  shows  how  the  modified  Cholesky  factorization  works  on  an  example  matrix  of 
size  n  =  4.  Consider  the  matrix 


0.3571 

-0.1030 

0.0274 

-0.0459 

-0.1030 

0.2525 

0.0736 

-0.3845 

0.0274 

0.0736 

0.2340 

-0.2878 

-0.0459 

-0.3845 

-0.2878 

0.5549 

The  eigenvalues  of  this  matrix  are  -.0767.  .1442  ,  .4004,  and  .9307  and  the  maximum  diagonal  element, 
y,  is  .5549.  Let  rj  =  r2  =  6.0555e  —  06,  which  is  the  value  of  machepsi  on  a  Sun  3/75,  using  double 
precision. 

In  the  first  iteration.  T44  is  the  maximum  diagonal  element,  therefore  row  and  column  4  are 
interchanged  with  row  and  column  1.  In  performing  the  test  of  whether  or  not  the  min  {.4„  —  -1  }  < 

;-ri<l<n  A” 

r-;.  the  minimum  occurs  at  .4;:  —  ^77  •  which  is  <  0.  and  consequently  the  algorithm  switches  to  phase 
2 

The  actual  lower  Gerschgorin  bounds  for  the  start  of  phase  2  are  [-.1633,  -.3086.  -.1548.  .1808], 
The  maximum  bound  is  the  bound  for  row  4.  hence  row  and  column  1  are  again  interchanged  with  row 
and  column  4,  and  because  this  bound  is  greater  than  0,  E\\  =  0. 

Prior  to  the  start  of  iteration  2,  the  updated  lower  Gerschgorin  bound  estimates  become 
[-.2564. -.1410. -.1401]  for  rows  2  through  4.  The  maximum  of  these  is  the  estimate  for  row  4.  re¬ 
sulting  in  a  diagonal  pivot  of  rows  and  columns  2  and  4.  The  actual  lower  Gerschgorin  bound  for  row  2 
is  -.1330.  therefore  £4 2  =  .1330. 

In  the  final  iteration,  the  eigenvalues  of  the  remaining  2x2  submatrix  are  .156329  and  -.052115. 
or  A i0  and  AAi,  respectively.  The  value  -A ,a  +  Aa,  -  A/ff)  is  .052119,  which  is  less  than  £22.  therefore 
the  algorithm  sets  both  £33  and  £44  to  the  value  added  to  A22  in  the  previous  iteration,  or  .1330. 

The  Cholesky  factors  and  pivot  vector  are 


0.5976 

0.0 

1 

-0.0769  0.8259 

0.1330 

4 

.  £  = 

and  P  = 

0.0458  -0.3442  0.4964 

0.1330 

3 

-0.1724  -0.4816  -0.1697  0.3082 

0.1330 

2 

where  PTAP  +  £  =  LLT .  and  P  is  1  permuted  by  the  transformations  recorded  in  P. 


The  r ;  1 1  i  o  hA'!)^/  —  A,  (.4).  where  Aj(.4)  is  the  most  negative  eigenvalue  of  A.  is  1.7.'!.  and  the 
condition  number  of  (.-1+  Ei  is  21. S.  In  comparison,  for  the  Gill,  Murray,  and  Wright  iTOSlj  algorithm. 

;  K  \-x_i  —  A>  (.4)  is  '.'.  IS.  and  the  condition  number  of  (.-1  -r  E)  is  39.2. 

4  Software  for  the  Modified  Cholesky  Factorization 

The  code  for  the  modified  Cholesky  factorization  is  a  straightforward  implementation  of  the  algorithm 
detailed  in  Appendix  l  of  Schnabel  and  Eskow  [1988].  It  is  organized  into  one  main  user-callable  sub¬ 
program  containing  three  smaller  subroutines,  each  of  which  are  called  only  once.  These  three  smaller 
subroutines  serve  to  initialize  variables  at  the  start  of  the  algorithm,  initialize  the  actual  Gerschgorin 
boutids  at  the  start  of  phase  2.  and  compute  the  factorization  of  the  final  2x2  submatrix  in  pha."-  2 
The  remainder  of  the  factorization  is  performed  by  the  main  subroutine.  In  particular,  while  the  rode 
for  pivoting  in  phase  1  and  phase  2  is  similar,  it  has  been  left  in-line  to  prevent  the  necessity  of  having 
Oin)  function  calls.  Because  the  row  and  column  pivoting  must  affect  only  th->  lower  triangle  of  the  input 
matrix,  this  code  is  lengthy  in  comparison  to  the  remainder  of  the  algorithm.  All  non-integer  variables 
in  the  code  are  double  precision. 

The  main  subprogram  is  called  by  modcholesky(ndun,n,A,g,macheps,ri,T2,pivot,E)  .  The  input 
parameters  to  this  subroutine  are: 

•  ndun  is  the  dimension  of  matrx  that  contains  A  in  the  calling  program. 

•  n  is  the  dimension  of  the  input  matrix  .-1. 

•  A  is  an  n  x  n  symmetric  matrix  (only  the  lower  triangular  portion  of  A,  including  the  diagonal,  is 
used,  and  it  is  overwritten  by  L). 

•  g  is  an  n  dimensional  work  vector. 

•  macheps  is  machine  epsilon. 

•  T\  is  the  reciprocal  of  the  tolerance  used  for  determining  when  to  switch  to  phase  2,  i.e.  1  /t\  is 
the  minimum  condition  number  of  a  positive  definite  input  matrix  which  may  be  perturbed  by  the 
algorithm. 

•  r2  is  the  tolerance  used  for  determining  the  maximum  condition  number  of  the  final  2x2  submatrix 
and  in  the  equation  for  Ejj- 
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TL-’  .nit put  parameters  arc 

•  L  is  -t  T' d  in  the  niT'.x  A  tin  the  lower  triangular  portion,  including  th<-  main  diagonn. 

•  pi’  ot  is  a  record  of  how  the  rows  and  columns  of  the  matrix  were  permuted  during  the  factorization 
That  is.  each  P,  is  initialized  to  i.  and  at  each  iteration,  if  rows  and  columns  i  and  j  are  switched, 
then  P ,  and  are  swapped. 

•  E  is  an  n-vector,  w  hose  i:h  element  is  the  amount  added  to  the  diagonal  of  .-1  at  the  i‘f'  iteration  of 
the  factorization. 

A  simpler  driver,  called  by  modchol(ndim.n.A .G.macheps.pivot.  E)  is  also  available.  This  driver  sets  the 
parameters  and  r-_  to  rnachf.ps'* .  and  the  remaining  input  and  output  parameters  are  identical  to  those 
for  modchoic >k’j. 

A  sample  drive-  program,  choldrtver.f.  and  its  output  are  included  with  the  code  T.e  modified 
Cho'esky  factcnzati  t).  The  driver  calls  a  function  macheps  to  compute  machine  epsilon  for  dcutde 
precision  arithmetic.  It  can  not  be  guaranteed  that  this  lunction  will  return  the  correct  value  of  niachcj' s 
cn  every  computer,  so  the  user  may  want  to  check  this  and.  if  necessary,  replace  the  call  to  macrif ps 
with  a  statement  that  assigns  the  actual  value  of  machine  epsilon  for  that  computer  to  eps .  The  driver 
program  also  calls  a  separate  subprogram,  mkmatnz.f.  to  generate  random  test  matrices  with  eigenvalues 
within  a  specified  range  Calls  to  both  rnodchoksky  and  modchoi  are  demonstrated  in  the  driver. 

Appendices  A  and  B  contain  the  sample  driver  and  sample  output,  respectively.  Appendix  C 
contains  the  modified  Cholesky  factorization  code 

Note  that  if  one  wishes  to  process  a  sparse  matrix  .1  incrementally  as  mentioned  in  Section  1. 
the  code  must  be  simplified  so  that-  all  pivoting  is  eliminated.  In  this  case  the  calculation  of  Gerschgorin 
bound  estimates  is  also  unnecessary  so  the  code  is  quite  simple.  The  diagonal  elements  must  still  be 
known  throughout  the  factorization,  but  the  rest  of  the  matrix  can  then  be  processed  incrementally,  with 
only  the  part  involved  in  the  current  elimination  step  needed  at  any  given  iteration. 


J 


References 


1'  Dennis.  J.  E  .  and  Schnabel,  R.  B.  Xumencal  Methods  for  Unconstrained  Optimization  and  .Son. 
lunar  Equations.  Preni . ir--— IIa.lI.  Englewood  Cliffs.  New  Jersey,  1983. 

21  Gill.  P  E.  and  Murray.  \V.  New;  n-type  methods  for  unconstrained  and  linearly  constrained  opti¬ 
mization.  Mathematical  Programming  28.  (1974).  311-350. 

[3]  Gill.  P.  E..  Murray.  W.  ano  Wright,  M.  II  Practical  Optimization.  Academic  Press,  London.  1981. 

_4;  Gould.  N  Private  communication.  1989. 

[5j  Schnabel,  R.  B.  and  Eskow.  E.  A  New  modified  Cholesky  factorization.  University  of  Colorado 
Department  of  Computer  Science  Technical  Report  Number  CU-CS-4 15-88.  (To  appear  in  SIAM 
J.  Sc i.  Stat.  Computing.) 


8 


o  o 


A  Sample  driver 


r  Din  it  for  run  in  t'iiti' d  chvli'-ky  factorization  algorithm, 

c 

integer  n.ndim 
double  precision  100,100) 
double  precision  Atwo(  100. 100) 
doxible  precision  g(1001 
double  precision  maxadd 
integer  pivot(1001 
double  precision  E(100) 
double  precision  eps.taul  ,tau2 
integer  7. 

double  precision  high. low 

ndim=  100 

i. '  machrps  subroutine  computes  machine  epsilon. 

1  the  following  line  may  be  replaced  by  assignment  to  eps 

(  '  of  correct  machine  epsilon  constant  for  yuur  machine 

call  machepsi.eps) 

C  Tolerances  used  by  modcholcsky  subroutine . 

C  tauJ  is  used  m  determining  when  to  switch  to  phase  2  and 

<'  tau2  is  used  in  determining  the  amount  to  add  to  the  diagonal 

C  of  the  final  2X2  submatnr. 

C  The  default  values  for  these  tolerances  can  be  used 

<’  by  calling  modchol  subroutine  instead  of  modcholcsky. 

C  The  default  values  for  taul  and  tau2  in  modchol  arc  :  eps  **  1/3. 

tau  1  =  eps  **  (../.l) 
tau'2  =:  eps  (1./3.) 

C  Initial  seed  for  random  .mber  generator  used  to  generate  lest 

C  matrices 

z  =  1000 

C  high  and  low  are  l  e  rany  of  the  eigenvalues  for  the  test  matrix 

C  to  be  generated. 

high  =  1.0 
low  =  -1.0 

The  first  test  problem  will  have  dimension  n—f,  so  that  the  entire 
problem  can  be  printed  out. 
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n  =  4 

print  ‘."TEST  PROBLEM  #1" 
print  ‘."Test  Matrix  of  size".n 

print  ‘."with  eigenvalues  within  the  range  of  ".low."  to  ".high 

call  mkmatrixf ndim.n.z.A .high, low, Atwo.g) 

print 

print  *, "Original  4X4  matrix" 
do  25  i=l,n 

25  print  (26).( A(i j) j  =  l,n) 

26  format  (4f20.8) 

call  modcholfndim.n.A .g.eps, pivot .E) 
print 

print  ‘."Matrix  after  factorization  with  1  m  the  lower  triangle" 
do  50  i=l.n 

■50  print  (26 ) .( A t i J ) j  =  l,n) 

print 

print  *.  "Iteration  Pivot  Amt  added  to  Aii" 
do  75  i=l,n 

75  print  ( 76 ) ,i .pivot ( i l.e(i) 

76  format  f i 2 . 1  Ox . i 2 . 1 0 x . f  1 2 .5 ) 

maxadd  =  E(nJ 
print 

print  ‘.".Maximum  amount  added  to  the  diagonal  is".ma.xadd 

C  The  next  3  test  problems  have  sixes  n— 35,50, &  75, 

C  with  eigenvalue  ranges  (~  1 ,  l].[- 1 , 1 0000],  &  [-10000,-1]  respectively. 

n  =  25 

print 

print  ‘/'TEST  PROBLEM  #2" 
print  ‘/'Test  Matrix  of  size",n 

print  ‘."with  eigenvalues  in  the  range  of  ",low,"  to  ".high 


call  mkmatrix(ndim,n,z, A, high, low, Atwo.g) 
call  modcho[(ndim,n.A,g,eps,pivot,E) 
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max  add  =  E(n) 


print  ’."Maximum  amount  added  to  th*’  diagonal  is".maxadd 


high  =  10000.0 
low  =  —1.0 
n  =  50 


call  mkmatrix(ndim,n.z.A.high.low,Atwo,g) 
print  *,"" 

print  ’."TEST  PROBLEM  #3" 
print  *,"Test  Matrix  of  size",n 

print  ’."with  eigenvalues  in  the  range  of  ".low,"  to  ".high 


call  modcholesky(ndim.n.A.g.eps.tau  1  ,tau2.pivot.E) 
maxadd  =  E(n) 

print  *, "Maximum  amount  added  to  the  diagonal  is", maxadd 

high  -  —1.0 
low  =  -10000.0 
n  =  To 

print 

print  ’."TEST  PROBLEM  #4" 
print  ’."Test  Matrix  of  size",n 

print  ’."with  eigenvalues  in  the  range  of  ".low."  to  ".high 
call  mkmatrix(ndim,n,z.A.high.low.Atwo.g) 
call  modchol(ndim.n,A,g,eps.pivot,E) 
maxadd  =  E(n' 

print  *, "Maximum  amount  added  to  the  diagonal  is", maxadd 


stop 

end 

£********#*********#*****************************#•##**#*#*####********** 
c  macheps 

c* ******************************************************  ***************** 
subroutine  macheps(eps) 

double  precision  eps 


130 


140 


macheps 


n 


double  precision  temp 
temp  =  1.0 

continue 

temp  —  temp  /  2.0 

if  ((1.0  +  temp)  .ne.  1.0)  goto  20 

eps  =  temp  *  2.0 

return 

end 
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B  Sample  Driver  Output 


TEST  PROBLEM  #1 

Test  Matrix  of  size  4 

with  eigenvalues  within  the  range  of  -1.0000000000000  to  1.0000000000000 


Original 

4X4  matrix 

0.35711021 

-0.10302945 

0.02737268 

-0.04594879 

-0.10302945 

0.25254612 

0.07358379 

-0.38451624 

0.02737268 

0.07358379 

0.23396662 

-0.28782367 

-0.04594879 

-0.38451624 

-0.28782367 

0.55494709 

Matrix  after  factorization  with  1  in  the  lower 

triangle 

0.59758699 

-0.10302945 

0.02737268 

-0.04594879 

-0.07689054 

0.82587804 

0.07358379 

-0.38451624 

0.04580534 

-0.34424172 

0.49639272 

-0.28782367 

-0. 17240912 

-0.48163633 

-0 . 16986202 

0.30827612 

Iteration 

Pivot  Amt 

added  to  Aii 

1 

1 

0 . 00000000 

2 

4 

0.13303961 

3 

3 

0.13303961 

4 

2 

0. 13303961 

Maximum  amount  added  to  the  diagonal  is  0.13303960618874 

TEST  PROBLEM  #2 

Test  Matrix  of  size  25 

with  eigenvalues  in  the  range  of  -1.0000000000000  to  1.0000000000000 

Maximum  amount  added  to  the  diagonal  is  1.2576119845957 


TEST  PROBLEM  #3 

Test  Matrix  of  size  50 

with  eigenvalues  in  the  range  of  -1.0000000000000  to  10000.0000000000 
Maximum  amount  added  to  the  diagonal  is  1.1271617927026 


TEST  PROBLEM  #4 

Test  Matrix  of  size  75 

with  eigenvalues  in  the  range  of  -10000.0000000000  to  -1.0000000000000 
Maximum  amount  added  to  the  diagonal  is  11618.452621394 
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C  Modified  Cholesky  Factorization  Code 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


subroutine  name  modcholesky 

authors  :  Elizabeth  Eskow  and  Robert  B.  Schnabel 
date  :  December,  19SS 

purpose  :  perform  a  modified  cholesky  factorization 

of  the  form  (Ptranspose)AP  +  E  =  L(Ltranspose), 
where  L  is  stored  in  the  lower  triangle  of  the 
original  matrix  A. 

The  factorization  has  2  phases: 

phase  1:  Pivot  on  the  maximum  diagonal  element. 

Check  that  the  normal  cholesky  update 
would  result  in  a  positive  diagonal 
at  the  current  iteration,  and 
if  so,  do  the  normal  cholesky  update, 
otherwise  switch  to  phase  2. 
phase  2:  Pivot  on  the  minimum  of  the  negatives 
of  the  lower  gerschgonn  bound 
estimates. 

Compute  the  amount  to  add  to  the 
pivot  element  and  add  this 
to  the  pivot  element. 

Do  the  cholesky  update. 

Update  the  estimates  of  the 
gerschgonn  bounds. 

—  largest  dimension  of  matrix  that  will  be  used 

—  dimension  of  matrix  .4 

—  n*n  symmetric  matrix  ( only  lower  triangular 
portion  of  A,  including  the  mam  diagonal,  is  used) 

—  n*l  work  array 
macheps  —  machine  epsilon 

iaul  —  tolerance  used  for  determining  when  to  switch  to 
phase  2 

tau2  —  tolerance  used  for  determining  the  maximum 
condition  number  of  the  final  2X2  submatnx. 

output  :  L  —  stored  in  the  matrix  A  (in  lower  triangular 


input  :  ndim 
n 

A 

9 


l&O 


200 


210 
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ooo 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c* 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


portion  of  A.  including  the  mam  diagonal I 

pivot  —  a  record  of  how  the  rows  and  columns 
of  the  malm  were  permuted  while 
performing  the  decomposition 

E  —  n*l  array,  the  ith  element  is  the 
amount  added  to  the  diagonal  of  A 
at  the  ith  iteration 

220 


mMwt*t******************t*******t*******^r*r*rt***^*rtt**r*mtt*»*er*rt 

subroutine  modcho)esky(ndim,n.A,gtmacheps,taul ,tau2.pivot.E)  modclioleskv 

integer  n.ndim 

double  precision  A(ndim.n),g(n),macheps,taul.tau2 
integer  pivot(n) 
double  precision  E(n) 


J 

lining 

imaxd 

i.itemp.jpl,k 

delta 

gamma 

norm] 

ming 

maid 

iaugamma 

phased 


—  current  iteration  number 

—  index  of  the  row  with  the  mtn.  of  the 

neg.  lower  Gersch.  bounds 

—  index  of  the  row  with  the  maximum  diag. 

element 

—  temporary  integer  variables 

—  amount  to  add  to  Ajj  at  the  jib  iteration 

—  the  maximum  diagonal  element  of  the  original 

matrix  A. 

—  the  1  norm  of  Afcolj).  rows  j-t-1 - >  n. 

—  the  minimum  of  the  neg.  lower  Gersch.  bounds 

—  the  maximum  diagonal  element 

—  iaul  *  gamma 

—  logical .  true  if  in  phased,  otherwise  false 


delta! .temp. jdmin.tdmm  —  temporary  double  precision  ears. 


230 


240 


integer  j.iming,i,imaxd,itemp jpl.k 
double  precision  delta, gamma 

double  precision  normj,  ming,maxd  :so 

double  precision  de!tal,temp jdmin,tdmin,taugamma 
logical  phase  1 

call  init(n,  ndim,  A,  phasel,  delta,  pivot,  g,  E, 
ming,  taul,  gamma,  taugamma) 

do  10  j  =  1,  n  — 1 

PHASE  1 

260 
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OOO  O  O  O  OOOu  O  O  O  ^  <-)  O  o 


if  (  phase  1  )  then 
C 

C  find  index  of  maximum  diagonal  dement  A(t.i)  where  i> —j 

C 

maxd  =  A  ( j  j ) 
imaxd  =  j 
do  20  i  =  j  +  1,  n 

if  (maxd  .It.  A(i.i))  then 
maxd  =  A(i,i) 

imaxd  =  i  zrz 

end  if 

20  continue 


pivot  to  the  top  the  row  and  column  with  the  max  dtag 
if  (imaxd  .ne.  j)  then 
swap  row  j  with  row  of  max  dtag 

260 

do  30  i  -  1 , j  —  1 
temp  =  A(j.i) 

A(j.i)  =  A(imaxd.i) 

A(imaxd.i)  =  temp 
0  continue 

swap  col)  and  row  maxdtag  between  j  and  maxdtag 

do  35  i  =  j  +  l,imaxd-l 

temp  =  A  ( i  j )  ;so 

A(ij)  =  A(imaxd.i) 

A(imaxd,i)  =  temp 
5  continue 

swap  column  j  with  column  of  max  dtag 

do  40  i  =  imaxd-fl,  n 
temp  =  A(ij) 

A(ij)  =  A  ( i ,  imaxd) 

A(i,  imaxd)  —  temp  300 

0  continue 

swap  dtag  elements 

temp  =  A(j  j) 

A(jJ)  =  A(imaxd,  imaxd) 

A(imaxd, imaxd)  =  temp 

swap  elements  of  the  pivot  vector 

310 
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Cj  u 


itcmp  =  pivot  (j) 
pivot(j)  =  pivot(imaxd) 
pivot  (imaxd)  =  itemp 

end  if 


C  Check  to  see  whether  the  normal  cholesky  update  for  this 

C  iteration  would  result  m  a  positive  diagonal. 

C  and  if  not  then  switch  to  phase  3'j: 

JPl  =  j+1 

if  (A(jj).gt.O)  then 

do  60  i  =  jpl,  n 

temp  -  ACi.j)  *  A(i.j)  /  A(j.j) 
tdmin  =  A(i.i)  —  temp 
if  (i  .ne.  jpl)  then 

jdmin  =  min(jdmin,  tdmin)  23.. 

else 

jdmin  —  tdmin 

end  if 

GO  continue 

if  (jdmin  .It.  taugamma)  phasel  =  .false. 

else 


phase  1  =  .false. 


end  if 


C 


if  (phasel)  then 

do  the  normal  cholesky  update  if  still  in  phase  1 


A(jvi)  =  dsqrt(A(jj)) 
do  70  i  =  j  +  1,  n 

A(ij)  =  A  ( i  j )  /  A  ( j  j ) 

“0  continue 

do  75  i=j+l,n 

do  80  k  =  j  +  1,  i 

A(i,k)  =  A(i,k)  -  (A(ij)  *  A(kj)) 
80  continue 

75  continue 

if  (j  .eq.  n— 1)  A(n,n)=dsqrt(A(n,n)) 


340 


350 


else 


360 
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c 

C  calcuh it<  llu  negatives  of  the  lower  gerschgorin  bounds 

( ' 

call  calcgersch! ndim.n.A.j.g) 

end  if 
end  if 


C 

C  PHASE  : 

C 

if  (.not.  phasel)  then 


C 

C 

c 


90 


if  (j  .ne.  n  — 1)  then 

find  (he  minimum  negative  gershgorin  bound 

do  90  i  =  j,n 

if  (i  ne.  j)  then 

if  (ming  gt.  g(i))  then 
ming  =  g(i) 
iming  =  i 
end  if 

else 

iming  =  j 
ming  =  g(j) 

end  if 
continue 


C 

c 

c 

c 

c 

c 

c 


pivot  to  the  top  the  row  and  column  with  the 
minimum  negative  gerschgorin  bound 

if  (iming  .ne.  j)  then 

swap  row  ]  with  row  of  min  gersch  bound 

do  100  i  =  1,  j  — 1 
temp  =  A(j,i) 

A  ( j ,  i )  =  A(iming,i) 
A(iming,i)  =  temp 
continue 

swap  col)  unlh  row  iming  from  j  to  iming 

do  105  i  =  j  +  l,iming-l 
temp  =  A(ij) 

A(i  j)  =  A(iming,i) 


350 


400 


410 
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O  O  O  O  O  O  OOO  O  O  O  Cj  o  o  o 


A(iming.i)  =  temp 
1  u5  continue 

C 

t'  >iiaf  rv/urnn  j  (nth  column  of  nun  gersch  bound 

C 

do  110  i  =  iming-H.  n 
temp  =  A(i  j) 

A(ij)  =  A(i.iming) 

A(i.iming)  =  temp 

110  continue  4t>o 

swap  diagonal  elements 

temp  =  A ( j  j ) 

A(jj)  =  A(iming.iming) 

A(iming.irmng)  =  temp 

swap  elements  of  the  pilot  vector 

itemp  =  pivot(j)  -i.i; 

pivot(j)  =  pivot(iming) 
pivot(inting)  =  itemp 

swap  elements  of  the  negative  gerschgonn  bounds  vector 

temp  =  g(j) 
g(j)  =  g(immg) 
g(iming)  =  temp 

end  if  mo 

Calculate  della  and  add  to  the  diagonal. 

delta  =  maz{0.  —  A(j,j)  +  mai{norm],laugamma}  ,delia_previous) 
where  normj-sum  of  \A(i.j)\.for  i=l,n. 

della^previous  is  the  delta  computed  at  the  previous  iteration, 
and  taugamma  is  taul*gamma. 


normj  =  0.0 

do  140  i  =  j+1,  n  450 

normj  =  normj  +  dabs(A(ij)) 

140  continue 

temp  =  max(normj,  taugamma) 
deltal  =  temp  —  A(jj) 
temp  =  0.0 

deltal  =  max(temp,  deltal) 
delta  =  max(deltal,  delta) 

E(j)  =  delta 

A(jj)  =  A(j  j)  +  E(j)  460 
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ooo 


c 

C  update  the  gerschgonn  bound  estimates 

C 

if  lA(j.j)  .no.  normj)  then 

temp  =  (normj/A(jj))  -  10 

do  150  i  =  j  +  1,  n 

g(i)  =  g(i)  +  dabs(  A(i.j))  *  temp 
150  continue 

end  if 

do  the  cholesky  update 

A  ( j  j )  =  dsqrt(A(j  j)) 
do  160  i  =  j  +  1 .  n 

A ( i j )  =  A ( i ,j '  /  A(j.j) 

continue 

do  165  l  =  j+1.  n 

do  170  k  =  j  +  1.  i 

A(i.k)  =  A(i.k)  -  ( A ( i J )  *  A ( k j ) ) 

continue 
continue 

else 

call  final2by2(ndim.  n.  A.  E.  j,  tau2,  delta, gamma) 
end  if 
end  if 

10  continue 

return 

end 

^■******************************************t*************************** 

C  subroutine  name  :  modchol 

C 

C  purpose  :  Simple  driver  for  the  modified  cholesky  algorithm, 

C  with  the  tolerances  set  to  the  default  values. 

C  t.e.  taul  =  tau2  =  macheps  **  1/3 

C 

C  input  :  n, ndim, A ,g, macheps 

C 

C  output  :  pivot, E 

C  (See  subroutine  modcholesky  above  for  details  on  all  parameters) 

Q****t$*******t*********t*******************tt******** ***************** 

subroutine  modchol( ndim, n, A, g, macheps, pivot, E) 


160 


170 

165 


500 


modchol 
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integer  ndim.  n 

double  precision  A(ndim.nbg(n).ma':h*’ps 
integer  pivot)  n) 
double  precision  E(n) 

double  precision  taul.tau2 
taul  =  macheps  **  (1./3.) 
tau2  =  taul 

call  modcholeskyfndim.n, A. g.  macheps. tau  1  ,tau2, pivot  F.) 


C* 

c 

c 

r 

i ' 
i ' 
i 

i ' 

( 

c 

c 

c 

c 

c 

c 


return 

end 


t  9  *  *  *  *  * 


subroutine  name  mit 

purpose  set  up  for  start  of  chol*  'ley  factorization 


input  .  n  u dun.  A.  taul 

output  phase!  —  boolean  value  set  to  true  if  in  phase  one. 
otherwise  false. 

delta  —  amount  to  add  to  Ajj  at  iteration  j 
pivot. g.  E  —  described  above  in  modchoksky 
nung  —  the  minimum  negative  gerschgortn  bound 

gamrna  —  the  maximum  diagonal  element  of  A 
taugamma  —  taul  *  gamma 


k*******9*********** 


subroutine  iiut(n.ndim. A, phas'd .delt a.pivot .g.E.ming. 
tau  1  .gamma. taugamma  I 


integer  n.ndim 

double  precision  A(ndim  m 

logical  phase  1 

double  precision  delta. g(n),E(n) 
integer  pivot(n) 

double  precision  ming.tau  1  .gamma, taugamma 


phasel  =  true, 
delta  =  0.0 
ming  =  0.0 
do  10  i=l,n 
pivot(i)  =  i 
g(‘)=  0 

E(i)  =  0 

10  continue 

c 


560 
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c  find  the  maximum  magnitude  of  the  diagonal  elements. 

.  if  any  diagonal  element  is  negative,  then  phasel  is  false. 

gamma  =  dabsf  A(  1 . 1 )) 

if  (Ai  1.1)  It.  0)  phasel  =  false. 

do  20  i  =  2.n 

if  (dabsf A( i.i) )  .gt.  gamma)  gamma=A(i.i) 
if  (A(i.i)  .It.  0)  phasel  =  .false. 

20  continue 

taugamma  =  taul  *  gamma 

c 

c  if  not  in  phasel.  then  calculate  the  initial  gerschgonn  bounds 

c  needed  for  the  start  of  phaseS. 

c 

if  (  . not. (phasel ))  call  calcgersch(ndim.n,A,l.g) 


5  70 


do  10  i  =  j,  n 
offrow  =  0.0 
do  20  k  =  j,  i  - 1 

20  offrow  =  offrow  +  dabsf  A(i.k)) 

do  30  k  =  i+1,  n 

30  offrow  =  offrow  +  dabs(A(k,i)) 

g(i)  =  offrow  —  A(i,i) 

10  continue 


610 
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('* 

c 

( ' 
c 
c 
c 
c 
c 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c* 


c 

c 


c 


c 

c 

c 


c 

c 

c 

c 

c 

c 


return 

end 


subroutine  name  :  finn!2by2 

purpose  Handles  final  2X2  subrnatnr  in  Phase  II. 

Finds  eigenvalues  of  final  2  by  2  submatnx. 
calculates  the  amount  to  add  to  the  diagonal, 
adds  to  the  final  2  diagonal  elements, 
and  does  the  final  update. 

input  :  ndim.  n.  .4.  E,  j.  tau2, 

delta  —  amount  added  to  the  diagonal  m  the 
previous  iteration 

output  :  ,-t  —  matnr  with  complete  l  factor  m  the  lower  tnanle. 

E  —  n*  1  vector  containing  the  amount  added  to  the  diagonal 
at  each  iteration. 

delta  —  amount  added  to  diagonal  elements  n—J  and  n. 


*******  *  *  *  *  *  • *  * *  *  *  *  ************************************************* 

subroutine  final2by2(ndim.  n.  A,  E,  j,  tau2.  delta. gamma) 

integer  ndim.  n.  j 

double  precision  Afndim.n),  E(n).  tau2.  delta.gamma 

double  precision  tl.  t2.  c3,lambdal  .lambda2.lambdahi.lambdalo 
double  precision  deltal,  temp 


find  eigenvalues  of  final  2  by  2  submatru 

tl  —  A( n—  1  ,n  - 1 )  t  A(ii.n) 
t2  =  A(n—  l.n—  1)  —  A(n,n) 
t3  =  dsqrt(t2*t2  +  •4.0*A(n,n-l)*A(n,n-l)) 
lambdal  =  (tl  —  t3)/2. 

Iambda2  =  (tl  +  t3)/2. 
lambdahi  =  max(lambdal,lambda2) 
lambdalo  =  min(lambdal,lambda2) 

find  delta  such  that: 

1.  the  12  condition  number  of  the  final 
2X2  submatnx  +  delta*!  <=  tau2 

2.  della  >=  previous  delta, 

3.  larnbdalo  +  delta  >=  tau2  *  gamma, 

where  lambdalo  is  the  smallest  eigenvalue  of  the  final 
2X2  submatnx 


final2by2 
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delta  1  =  !  larnbdahi  -  lambdalo)/(  1 .0  —  tau21 
deltal=  ma.\(deltal  .gamma) 
del'. al=  tau2  *  deltal  —  lambdalo 
temp  =  0.0 

delta  =  max(delta.  temp) 
delta  =  max(deltal.  delta) 

if  (delta  gt.  0.0)  then 

A(n—  1  ,n  —  1 )  =  A(n—  1  ,n  —  1 )  +  delta 
A(n,n)  =  A(n,n)  +  delta 
E(n— 1)  =  delta 
E(n)  =  delta 
end  if 

final  update 

A(n—  I.n  —  1)  =  dsqrt(  A(n— 1  ,n— 1 )) 
A(n.n-l)  =  Afn.n  — l)/A(n-l.n—  1) 

A(n.n)  =  A(n.n)  -  ( A(n.n—  l)*A(n.n—  1)) 
A(n.n)  =  dsqrt(A(n.n)) 

return 

end 
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